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ABSTRACT 

We find general relativistic solutions of equilibrium magnetic field configurations in 
magnetars, extending previous results of Colaiuda et al. (2008). Our method is based 
on the solution of the relativistic Grad-Shafranov equation, to which Maxwell's equa- 
tions can be reduced in some limit. We obtain equilibrium solutions with the toroidal 
magnetic field component confined into a finite region inside the star, and the poloidal 
component extending to the exterior. These so-called twisted-torus configurations have 
been found to be the final outcome of dynamical simulations in the framework of New- 
tonian gravity, and appear to be more stable than other configurations. The solutions 
include higher order multipoles, which are coupled to the dominant dipolar field. We 
use arguments of minimal energy to constrain the ratio of the toroidal to the poloidal 
field. 
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1 INTRODUCTION 

In this paper we construct models of non rotating, strongly 
magnetized neutron stars, or magnetars (Duncan & Tliomp- 
son 1992), in general relativity. We extend our previous work 
(Colaiuda et al. (2008)) based on a formalism developed in 
Konno et al. (1999) and in loka & Sasaki (2004), includ- 
ing the toroidal field in a twisted-torus configuration. The 
extension to this field geometry is accomplished with an ap- 
propriate choice of the function which determines, point by 
point, the ratio between the toroidal and poloidal compo- 
nents of the magnetic field. The non-linear relation among 
the functions defining the toroidal and poloidal fields, nat- 
urally leads to couplings between different multipoles, thus 
making inadequate the one multipole solution which is usu- 
ally assumed. 

A motivation to find consistent equilibrium solutions 
in general relativity with this particular geometry comes 
from recent progresses in numerical magneto-hydrodynamic 
(MHD) simulations, that have made possible to study the 
dynamics and the stability of magnetic stars. By following 
the time evolution of generic initial configurations, in the 
framework of Newtonian gravity and using polytropic equa- 
tions of state, Braithwaite & Spruit (2004) and Braithwaite 
& Nordlund (2006) (see also Braithwaite & Spruit (2006)) 
have found magnetic field configurations, which are stable 
on timescales much longer than the Alfven time: they de- 
cay only due to finite resistivity. These configurations are 
roughly axisymmetric; the poloidal field extends throughout 
the entire star and to the exterior, while the toroidal field 



is confined in a torus-shaped region inside the star, where 
the field lines are closed. These configurations were named 
twisted-torus. Furthermore, Yoshida et al. (2006) have shown 
that such configurations are not significantly affected by ro- 
tation; Geppert & Rheinhardt (2006) studied the depen- 
dence of magnetostatic equilibrium configurations on the ro- 
tational velocity and on the initial angle between rotation 
and magnetic axis, finding hints for the existence of a unique 
stable dipolar magnetostatic configuration, independent of 
the initial field geometry. 

We must remark that this particular field geometry re- 
sulting from dynamical simulations is obtained assuming 
that outside the star there is vacuum; consequently, out- 
side the star electric currents are forbidden and the magnetic 
field can only be poloidal. This implies that the toroidal field 
cannot extend to the exterior and that the field lines which 
cross the surface are purely poloidal, whereas the field lines 
confined inside the star can maintain a mixed (poloidal and 
toroidal) structure. This configurations appear to be stable 
on dynamical timescales, probably due to magnetic helic- 
ity conservation, which requires the persistence of a toroidal 
component of the field. Notice, however, that different so- 
lutions including a magnetosphere may be possible. In this 
case, the toroidal field could also extend to the external re- 
gion leading to a twisted magnetosphere (Lyutikov 2006; 
Pavan et al. 2009). 

In our perturbative approach, there is a free parame- 
ter which represents the ratio between the toroidal and the 
poloidal components of the magnetic field. We estimate the 
value of this parameter, by identifying the configuration of 
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minimal energy at fixed magnetic helicity. We also mention 
that in our configurations the external field has mainly a 
dipole structure, with small corrections from higher multi- 
poles, and that, confirming a previous suggestion by Pren- 
dergast (1956), the toroidal and poloidal fields have ampli- 
tudes of the same order of magnitude. Similar configurations 
have been found in Newtonian models including rotation, in 
Yoshida & Eriguchi (2006), Yoshida et al. (2006). 

In this paper we consider non-rotating stars because 
observed magnetars have a very slow rotation rate, although 
high rotation rates may occur in the early stages of their 
evolution. 

The structure of the paper is the following. The model 
is presented in Section 2. In Section 3 we discuss a config- 
uration with purely dipolar magnetic field, neglecting the 
couplings with higher multipolcs. In Section 4 we include 
the I = 1 and / = 2 field components; Section 5 accounts 
for the general case, including all multipoles and their cou- 
plings. In Section 6 we compute the total energy and the 
magnetic helicity, and estimate the parameter C,^ which de- 
termines the ratio between toroidal and poloidal fields by 
energy minimization; we also compute the magnetic energy, 
and compare the contributions of the poloidal and toroidal 
fields for different values of i^o. In Section 7 we discuss the 
results and draw the conclusions. 



2 BASIC EQUATIONS AND FORMALISM 

We Eissume that the (non-rotating) magnetized star is sta- 
tionary and axisymmetric. We further assume that the mag- 
netic field acts as a perturbation of a spherically symmetric 
background describing a spherical star. The magnetized fluid 
is described within the framework of ideal MHD, in which 
the effects of finite electrical conductivity axe neglected. Rig- 
orously speaking, this approximation is only valid while the 
crust is still completely liquid and while the core matter 
has not yet performed the phase transition to the superfluid 
state, which is expected to occur at most a few hours after 
birth (see e.g. section 5.1 in (Aguilera et al. 2008) and ref- 
erences therein). The onset of superfluidity and/or crystal- 
lization limits the period during which magnetostatic equi- 
librium can be established. Both the melting temperature 
and the critical temperature of transition to the superfluid 
state, are between 10® to 10^^° K, and a typical neutron 
star quickly cools down below this temperature in a few 
hours. However, since the characteristic Alfvon time is on 
the order of ta « 0.01 — 10 s, depending on the background 
field strength, there is ample time for the magnetized fluid 
to reach a stable state, as shown in Braithwaite & Spruit 
(2006), while the state of matter is still liquid. After the 
crust is formed, the magnetic field is frozen in, and it only 
evolves on a much longer timescale due to Ohmic dissipa- 
tion or, in some case, due to the Hall drift (Pons & Geppert 
2007). Therefore, it is reasonable to expect that the MHD 
equilibrium configurations set within the first day after for- 
mation, will fix the magnetic field geometry for a long time. 

Here we first summarize the basic equations of ideal 
MHD in the framework of general relativity and then intro- 
duce the perturbative approach. Next, wo obtain the form 
of the electromagnetic potential in the case of twisted-torus 
configurations, and derive the relativistic Grad-Shafranov 



equation. In what follows no assumption is made on merid- 
ional circulation; thus, in principle, the meridional flow 
may be non-vanishing. We use spherical coordinates, x'^ = 
{t,x°',4>), where a;" = {r,d). A stationary axisymmetric 
space-time admits two killing vectors, r) = d/dt and ^ = 
d/d<j}, and with our coordinate choice all quantities (includ- 
ing the components of the metric tensor g^^) are indepen- 
dent of t and d>. 



2.1 Equations of ideal MHD in General Relativity 

According to a comoving observer with four- velocity w*^ , the 
stress-energy tensor of a perfect fluid with an electromag- 
netic fleld is 



where 



■ -'■ fluid "r -'em > 



^2 _ gi,gu 



(1) 

(2) 

(3) 



As usual, Euler's equations are found by projecting the equa- 
tion Ti^"^^ = orthogonally to u^^ 



{p + P)ai, + -I- Ui,uP,. -U=0 



(4) 



-eaf3j5 ■ 



(5) 



where /^j = F^vJ^ is the Lorentz force and a^j = u'^u^-,^ . 
Here, F^^ = duAp, — dfiA^ is the Mcixwell tensor, in terms 
of which the electric and magnetic fields can be defined as 

1 

2' 

The basic equations of ideal, general relativistic MHD 
are, then: (i) the continuity equation [nu'^)-.^ = 0; (ii) 
Maxwell's equations F'^",, = AnJ'^; (iii) the condition of a 
vanishing electric field in the comoving frame Ejj, = F^^u" = 
and (iv) Euler's equations (4). 

2.2 The perturbative approach and the form of 
the electromagnetic potential 

We treat the magnetic field as an axisymmetric perturbation 
of a spherically symmetric background and seek for station- 
ary solutions. The background metric is 

ds^ = -e'-Mdt^ + e^Wrfr^ + r\de^ + sin^ 0dc^>^') , (6) 

where !/(r), fi{r) arc solution of the unperturbed Einstein 
equations (the TOV equations) for assigned equations of 
state. The unperturbed 4-velocity is m'' = (e""/^, 0, 0, 0). To 
model the unperturbed neutron stax we use the equation of 
state of Akmal, Pandharipande and Ravenhall called APR2 
(Akmal at al. 1998), with a standard equation of state for 
the stellar crust (see Benhar et al. (2004)), which results in 
a neutron star of mass M = 1.4 and radius R = 11.58 
km. 

It can be shown (sec for instance Colaiuda et al. 
(2008)) that (F^^, A^,, J'') are of order 0(B), and the per- 
turbations {6u'^, 5p, 6P, 5n, Sg^^, SGfi,^, ST^^) are of order 
0{B^). Therefore, at first order in B the magnetic fleld 
is coupled only to the unperturbed background metric (6), 
whereas the deformation of the stellar structure induced by 
the magnetic field, which we do not consider in this paper, 
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appears at order O(B^). Furthermore, (S*, At, J*, Fti,) = 
0{B'') and (/,, U) = 0{B*). 

With these assumptions, the potential at 0{B), has 
the form A^{r,d) = {0,Ar,Ao,A^). With an appropriate 
gauge choice we can impose Ag = and write the potential 
as 



yl^ = (0,e— S,0,V) , 



(7) 



where the components of A^ are expressed in terms of two 
unknown functions, E(r, 6') and xp{r,0). 

A further simplification of A^ is possible by exploit- 
ing the fact that U = -^,rJ'' - tp,eJ^ = 0{B^). Using 
Maxwell's equations and neglecting higher order terms, we 
find 



(8) 



where ip = sin0 S_e. This result implies V' = V'(V') and allows 
us to write 



sin6»S,e = CWV' , 



(9) 



where f (V') is a function of tp of order 0(1). 

The function C, represents the ratio between the toroidal 
and poloidal components of the magnetic field; difi^erent 
choices of this function lead to qualitatively different field 
configurations. The simplest case is C = constant (loka & 
Sasaki 2004; Colaiuda et al. 2008; Haskell et al. 2008), but 
with this choice (like with other simple choices of Ci^-')) the 
field lines of the toroidal field reach the exterior of the star, 
where there is vacuum. However, the magnetic field in vac- 
uum can only be poloidal (see, for instance, Colaiuda et 
al. (2008)), thus this solution presents an inconsistency. To 
avoid this inconsistency, one should consider a non- vacuum 
exterior, i.e. a magnetosphere, but modelling a neutron star 
magnetosphere is a quite difficult task, especially in gen- 
eral relativity. An alternative choice is to assume that the 
magnetic field is entirely confined inside the star (loka & 
Sasaki 2004; Haskell et al. 2008), but in this way the param- 
eter ( must assume particular, ad-hoc values; or, one can 
instead accept that the toroidal field has a discontinuity at 
the stellar surface, vanishing in the exterior (Colaiuda et al. 
2008); in this way the entire range of ^ can be studied, but 
the discontinuity in the field will, for consistence, imply the 
existence of ad-hoc surface currents. 

A different choice is made in this paper, where we as- 
sume that the toroidal field is confined in a toroidal region 
inside the neutron star, such that its field lines never cross 
the stellar surface, as in the twisted-torus configuration. As 
mentioned in Section 1, Newtonian numerical simulations 
(Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006; 
Braithwaite & Spruit 2006) suggest that these configura- 
tions are indeed a quite generic outcome of the evolution of 
strongly magnetized stars. 

The twisted-torus configuration can be obtained by 
choosing the following form of the function f 



= Co - 1] edv/v^l - 1) 



(10) 



A similar choice has been made, in a Newtonian framework, 
in Yoshida et al. (2006). In equation (10), is a constant of 
order 0(1); is a constant of order 0{B): it is the value of 
xjj at the boundary of the toroidal region where the toroidal 
field is confined (this boundary is tangent to the stellar sur- 
face); finally, Q{\tp/'ip\ — 1) is the usual Heaviside function. 



With this choice, the function vanishes at the stellar sur- 
face, where r = R, and the magnetic field 



r2 sin-^ * 



ed^/V'l-i) (11) 



has no discontinuities. 



2.3 The relativistic Grad-Shafranov equation 

The Grad-Shafranov (GS) equation, which allows to deter- 
mine the magnetic field configuration, can be derived from 
the 0-component of Maxwell's equations 



1 



4-jYf' 



[ip,ee - cot etp,e] (12) 



and from the a-components of Euler's equations (4), as fol- 
lows. Euler's equations give 

fa = (P + P)aa + P,a + UaU" P,^ 

= {p + P){^-e^5u*) ^ + P,a + 0{B^) . (13) 

For barotropic equations of state P = P{p), the first princi- 
ple of thermodynamics allows to write 



Pa = {p + P)(ln 
then (13) yields 



p + P 



fa = {p + P)x,< 



(14) 



(15) 



where x = x(^)^)- On the other hand, the o-components 
of the Lorentz force = F^^.^J'^ can be written as (see 
Colaiuda et al. (2008)) 



fa = - 
r 

where, in the present case. 



(16) 



Therefore, 



'_ ^0 

47r 



^%-J^ = (p + P)X,a. 

r*" sin 6/ 

From x,re — Xfir = it follows that 



- i',e 



(p + P)r^ sin^ e 



(17) 



, 



■th 

' (p + P)r2 sin^ 
which implies 

( ^ 

\ip + P)r^ sin^ e 

with Co, ci constants of order 0{B), 0(1) respectively. 
Hence, turns out to be 



^ = F(V') = co + ciV -h O(B^) , (18) 



= ^^[V'-3viv/v5| + 2v^vvi'] e(iv/v;i-i) 

+{p + Py sin^ e [co -I- Clip] . 



(19) 
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From Eqns. (12), (19) the rclativistic GS equation at first 
order in B is finally obtained; 



47r 



ip H — ip 



_e Co 

47r 



1.3 ,T2 



[V'.ee - cot^V.e] 



= {p + Py sin e [co + ci^P] . (20) 

If we now define i/)(r, 6) = sin da{r, 0),e and expand the func- 
tion a(r,9) in Legendre polynomials 



a(r,e) = '^ai(r)Pi{cose) 



(21) 



the GS equation rewrites as 



sm» / _A // , 



z/'-A' , Z(Z + 1) 



ai 



3 



e 



-J aiPie sind 



1 d 



y~^aiPi,flSing 



OiPi,e sin6|o('P(',9 sin( 



+ I 0(Oi/ai//P(,eP!',9Pz",9sin^6l 



: (p + P)r^ sin^ 61 



Co + ci ^OiP(,e sin^ 



(22) 



Finally, projecting Eq. (22) onto the different harmonic com- 
ponents, we obtain a system of coupled ordinary differen- 
tial equations for the functions 0((r). The projection is per- 
formed using the property 



21' + 1 



Pi, g Pi,, g sin 9 (le 



(23) 



+ 1) Jo 

If we consider the contribution of n different harmonics, we 
need to solve a system of n coupled equations, obtained from 
(22), for the n functions ai{r). 

2.4 Boundary conditions 

The functions az(r) must have a regular behaviour at the 
origin; by taking the limit r — > of the GS equation one can 
find 



ai{r ^ 0) = air 



(24) 



where ai are arbitrary constants to be fixed. 

Outside the star, where there is vacuum and the field is 
purely poloidal, equations (22) decouple, and can be solved 
analytically for each value of I. The solution can be ex- 
pressed in terms of the generalized hypcrgcornctric functions 
(P([Zi, ^2], [I3], z)), also known as Barnes' extended hyperge- 
ometric functions, as follows: 

ai = Air-^ F{[l,l + 2],[2 + 2l],z) 

+A2r'+'' F{[l-l,-l-l],[-2l],z) . (25) 

where z — 2M/r and Ai and A2 are arbitrary integration 
constants, which must be fixed according to the values of 
the magnetic multipole moments. Regularity of the external 



solution at r = 00 requires A2 = for all multipoles. For 
example, for I = 1, 2, 3 we have 



ai oc 



02 oc 



as oc 



ln(l -z) + z+' 



(4 - 3z) ln(l -z)+4:z-z'^-^ 



(15 - 2O2 + 6«^) ln(l - z) + 15z 



2 ^ ^ 12 



(26) 



At the stellar surface we require the field to be continuous. 
This condition is satisfied if ai and a'l are continuous. For 
practical purposes, the boundary conditions at r = P can 
be written as 



at 



'r 



(27) 



where /( is a rclativistic factor which only depends on the 
star compactness 2M/R (in the Newtonian limit all /; = 
1), and can be numerically evaluated with the help of any 
algebraic manipulator. For our model {2M/R = 0.357), the 
values of /; for the first five multipoles are 1.338, 1.339, 
1.315, 1.301, and 1.292, respectively. 

In general, there are n + 2 arbitrary constants to be 
fixed: the n constants ai, associated to the condition (24), 
Co and ci . Thus, wo need to impose n+2 constraints, of which 
n + 1 are determined by the boundary conditions, n condi- 
tions are provided by Eq. (27), i.e. by imposing continuity in 
r = P of the ratios a'l/ai. The overall normalization of the 
field gives another condition, which is fixed by imposing that 
the value of the I = 1 contribution at the pole is Bpde = 10^^ 
G (this corresponds to set ai(P) = 1.93- 10~^ km). The rea- 
son for this choice is that the surface value of the magnetic 
field is usually inferred from observations by applying the 
spin-down formula, and assuming a purely dipolar external 
field; for magnetars, the value of B estimated in this way is 



10^ 



10 G. The remaining condition will be imposed 



as follows. 

In the case of a purely dipolar field {I = 1), we shall 
assume ci — 0. In the general multipolar case, we choose 
to impose that the external contribution of all the I > 1 
harmonics, i.e. X)(>i '^i(-R)^' is minimum. 



3 THE CASE OF PURELY DIPOLAR FIELD 

We begin discussing the simplest case of a purely dipolar 
configuration, in which all couplings with higher order mul- 
tipoles are neglected in Eq. (22) (oi>i = 0). In this case, for 
any assigned value of Co there exists an infinite set of solu- 
tions, each corresponding to a value of ci; these solutions 
describe qualitatively similar magnetic field configurations. 

However, when higher order harmonics are taken into 
account, as we will see in the next Section, the picture 
changes. For instance, when Co = and the I = 1,2 har- 
monic components are included, the equations for ai and 02 
decouple: the equation for oi is the same as in the purely 
dipolar case, but a solution for 02 satisfying the appropri- 
ate boundary conditions exists only for a unique value of ci . 
Therefore, in the general case ci is not a truly free parameter 
(this is true also for Co 0), and the fact that in the purely 
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Figure 1. The profiles of the tetrad components of the magnetic field = 0), B^g^{6 = 7r/2), B(^^^{6 = 7r/2)) are shown for the 

purely dipolar case with Co = km~^, Co = 0.40 km~^ and Co = 0.80 km~^. 



dipolar case it looks as such, is an artifact of the truncation 
of the I > 1 multipoles. In order to provide a mathemati- 
cally simple example, which will be useful to understand the 
structure of the twisted-torus configurations, in this Section 
we shall consider the simplest case Ci = 0. 

By projecting Eq. (22) onto the I = 1 harmonic, and 
neglecting all contributions from I > 1 terms we find 

1 ( -X „ ^ -x ^' , 2 



--^1 (3/4)9 



-ai sin 



- 1 



ai + 3ai 



-ai sm 



2a? sin"* 9/^^ 



= (3/4) / co(p + Py sin^ Ode = co(p + P)r^ . (28) 
Jo 

The tetrad components of the magnetic field (i.e. the com- 
ponents measured by a locally inertial observer) are: 



0.5 



0.0 




Figure 2. The projection of the field lines in the meridional plane 
are shown for the purely dipolar case with Co = 0.40 km~^. The 
toroidal field is confined within the marked region. 



B, 



(r) 



sin 6 

_ A 
e 2 



B 



r sm 

e 



■i',r 



(4>) = 



*CoV' (l^/V-l-i) 



• edv^/VJl - 1) , 



(29) 



where = — ai sin^ 6. 

The profiles of the tetrad components of the field inside 
the star, are plotted in Fig. 1 for increasing values of ^o; 
B(r) is evaluated in {6 — 0) and B^e^, 5(0) in {6 = 7r/2). In 
Fig. 2 wo show the projection of the field lines in the merid- 
ional plane, for C,o = 0.40 km~^. Figs. 1 and 2 show that the 
toroidal field 5(0) is confined within a torus-shaped region; 
its amplitude ranges from zero, at the border of the region, 
to a maximum, close to its center. At the stellar surface and 
in the exterior B'^ vanishes, showing that there is no dis- 
continuity in the toroidal field. The panels of Fig. 1 show 
the field profiles for different values of Co: larger values of 
Co correspond to a toroidal field with increasing amplitude, 
confined in an increasingly narrow region close to the stel- 
lar surface, while the amplitude of the poloidal components 
(i?(^), B(e)) decreases. We remark that this implies that in- 
side the star we cannot have a twisted-torus configuration 
where the toroidal component dominates with respect to the 
poloidal one: if |-B(,^)| becomes larger with respect to |-B(r)| 
and |S(e)|, the domain where it is non vanishing shrinks. 



4 THE CASE WITH / = 1 AND I = 2 
MULTIPOLES 

Wc now proceed with our investigation considering the 1 = 1 
and / = 2 contributions, and setting o;>2 = 0. The projec- 
tion of the GS equation (22) onto the harmonics / = 1 and 
I = 2 gives the following coupled equations 



1 / „ , ^z.'-A'_, 2 



47r 



e ai + e 



47r 



/^3/4)e( 



-ai joi 

2 



-ai — 3a2 cos 6 . 2 
— sm ( 



Co 



— oi — 3o2 cos ^ -h 3(oi -h 3o2 cos 6) 

+ 2 sin* ^ ( — oj — 9oi02 cos ( 



-ai — 3a2 cos 9 . 2 
— sm 



-27oiai cos^ e - 27al cos^ 9) jil? 



{p + P)r^ ( CO - ^cioi ) , 



(30) 
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Figure 3. The projection of the field lines in the meridional plane are shown for = km ^ and 02 (-R) /ai (-R) = 1, 1/2, 1/4 respectively, 
and for a(>2 = 0- The dashed line corresponds to V = 0- 



, -A " , — 



A' , 6 

02 — -^0L2 



-ai — 3a2 cos 9 . 2 
— sm ( 



ai — 3a2 cos 6 + 3(ai + 3o2 cos ff) 



— ai — 3a2 cos^ . 2 „ 
— sm 



+ 2 sin'' 6( — o? — 9o?02 cos ^ 



-27aiai cos^ 6 - 27al cos^ 6*) /^^ 



~(p + P)r'^cia2 



(-3cos6'sm^6') 



(31) 



We integrate tliis system by imposing the boundary condi- 
tions discussed above, i.e. a regular behaviour at the origin 
(equation (24)) and continuity at the surface of ai, a[, a2, 02 
with the analytical external solutions given by (26). 

Let us first consider the simple case Co = 0. Eqns. (30), 
(31) decouple, and become 

A' 



-X II , — I 2 

e Oi+e — - — ai -ai 

2 r-^ 



Co 



:C\ai 



—\ II I 
e 02 + e 



. I/' - A' 



-a2 



-02 



167r , , 2 
- — (p+p)r C102 . 



(32) 



There are four constants to fix (qi, a2, co, ci) and three con- 
ditions: ai{R) — 1.93 ■ 10^"^ km (normalization) and the ra- 
tios a'i{R)/ai{R) and a'2{R)/a2{R) from the matching with 
the exterior solutions; thus, we need an additional require- 
ment. We remark that we cannot impose ci = as in the 
purely dipolar case, because the ratio a^i^R) / a2[R) depends 
only on ci, and the matching with the exterior solution is 
possible only for a particular value of ci, i.e. ci = 0.84 km~^. 

If we impose that |o2(i?)| is minimum, we find that this 
condition yields the trivial solution 02 (r) = (with non- 
vanishing ai). Indeed, from Eqns. (32) it is straightforward 
to see that 02(7-) = is a solution of the system. When 
Co 0, equations (30), (31) are coupled, but they still allow 
the trivial solution 02(7") = 0, which minimizes |o2(J?)|, with 
non- vanishing ai . The existence of this solution is a remark- 
able property of this system, and it is due to the fact that 
the integral in 6 on the left-hand side of Eq. (31) vanishes 



for 02 = (the integrand becomes odd for parity transfor- 
mations ^ — » TT — 6). Hence, if we look for a solution which 
minimizes the contributions from the Z > 1 components at 
the stellar surface, we have to choose the trivial solution 
a2(r) = 0. 

If, instead, we do not require that a2{R) is minimum, 
and we assign a finite value to the ratio 02 (-R) /ai (7?) , we find 
a non-trivial field configuration which is non symmetric with 
respect to the equatorial plane. This feature is shown in Fig. 
3, where the projection of the field lines in the meridional 
plane are plotted for Co = and a2{R)/ai{R) equal to 1, 
1/2 and 1/4, respectively. 



5 THE GENERAL CASE 

When all harmonics are taken into account, there exist two 
distinct classes of solutions: those symmetric (with respect to 
the equatorial plane) , with vanishing even order components 
{a2i = 0), and the antisymmetric solutions, with vanishing 
odd order components (02;+! = 0). Both solutions satisfy 
the GS equation (22). Let us consider the symmetric class. 
When 02i = the integrals arising when equation (22) is 
projected onto the even harmonics, which couple odd and 
even terms, vanish since the integrands change sign under 
parity transformations. Therefore, the symmetric solutions 
can be found by setting 02; = 0, projecting Eq. (22) onto the 
odd harmonics and solving the resulting equations for 02i+i. 
Similarly, the integrals in equation (22) projected onto the 
odd harmonics, vanish when a2i+i — 0; thus, we can consis- 
tently set a2i+i = 0, and find the antisymmetric solutions 
using the same procedure. 

In Section 4, we set the value of ai at the surface to be 
1.93-10"^ km and we minimized the 1 = 2 contribution. It is 
clear that, since the I = 1 and I = 2 iiiultipolcs belong to dif- 
ferent families, any attempt to minimize the relative weight 
of one with respect to the other leads to the trivial solution. 
The properties of equation (22) discussed above, tell us that 
if ai (R) 7^ we cannot consistently set to zero the remaining 
odd order components a2i+i- However, we have the freedom 
of setting to zero all even terms 02; . Therefore, since we have 
chosen to minimize the contributions of the I > 1 harmon- 
ics outside the star, we shall focus on the symmetric family 
of solutions {a2i = 0); we will briefly discuss an example 
belonging to the antisymmetric family in subsection 5.4. 
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Figure 4. The profiles of the tetrad components of the magnetic field = 0), B(^g^(6 = 7r/2), B(^^j(6 = 7r/2)) are shown for Co = 

km-i, Co = 0.40 km-i and Co = 0.80 km"!, and 1 = 1,3. 





Figure 6. The projection of the field lines in the meridional plane is shown for Co = km^^ and I = 1,3. The left panel refers to the 
solution corresponding to the absolute minimum of |a3(-R)/ai(-R)|; in this solution i/i has no nodes. The center and right panels refer 
to solutions corresponding to relative minima of \a3{R)/ai{R)\; in these cases tp has one and two nodes, respectively. The dashed lines 
corresponds to the ip = surfaces. 
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Figure 7. The profiles of the tetrad components of the magnetic field (B(r) = 0), {6 = tt/S), B(<^) {6 = 7r/2)) for the case including 
; = 1, 3, 5, with Co = km-i, Co = 0.61 km"! and Co = 3.00 km"!. 
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Figure 8. The projection of tlie field lines in the meridional plane is shown for Co = km 
respectively, and for i = 1, 3, 5. The toroidal field is confined within the marked region. 



0.61 km-i and Co = 3.00 km" 



5.1 The case with I = 1 and I = 3 

We now consider the system of equations including only the 
/ = 1 and / = 3 components. The projected system is 



A' , 2 

Oi jOl 



/ (3/4) Co - SVIV'/V'I + 2V'VV;' 
•e(l'(/-'/'0| - I)sin6ld6l 



Co — -ci I ai — -as 



(33) 



of the truncation in the harmonic expansion, and disappears 
as higher order harmonics are included. 

For completeness we also mention that the solutions cor- 
responding to the local minima of \a'j{R) \ different from the 
absolute minimum, correspond to very peculiar field con- 
figurations (see Fig. 6). The function ip has nodes on the 
equatorial plane, therefore the field lines lie in disconnected 
regions; for a fixed value of ^o, the number of nodes increases 
as I as (J?) I increases. These peculiar solutions exist for any 
value of Co, and appear also when higher order harmonic 
components are considered. Thus, they are not artifacts of 
the /-truncation. 



1 / -A ,/ , -x i^'-y , 12 

+ ^ ^"(7/48) Co (V- - 3V|V'/V'I + 2^^^' 
■Q{\ip/'4>\ - 1)(3- 15cos^ e)amede 

lo 



where 



-ai + 



a3(3- 15cos^ 61) 



(34) 



(35) 



We again impose regularity at the origin (Eq. (24)), conti- 
nuity in r = _R of ai,a'i,a-j,a'[i with the vacuum solutions 
for ai(r), a3(r) given by Eq. (26), and we fix ai{R) = 
1.93 • 10~^ km by normalization. For the remaining con- 
straint we choose the solution that minimizes the absolute 
value of a3{R). We find that there is a discrete series of local 
minima of \a3{R)\, and we select among them the absolute 
minimum. 

Fig. 4 shows the profiles of the tetrad field components 
(see Eq. (29)) obtained by numerically integrating Eqns. 
(33), (34), for different values of Co- S(r) is evaluated at 
6 = 0, while -6(9), -8(0) arc evaluated at ^ = tv/2. As ("o in- 
creases, the magnitude of the toroidal field becomes larger, 
but the region where it is confined shrinks, as already found 
in section 3. The projection of the field lines in the merid- 
ional plane is shown in Fig. 5 for the same values of (^q. It 
shows that, for ^ 0.40 km~^, the magnetic field lines lie 
in disconnected regions, separated by dashed lines in the fig- 
ure. Inside these regions, the function tp has opposite sign 
and no toroidal field is present. A similar phenomenon has 
been discussed in Colaiuda et al. (2008). As we will see in the 
next Section, the occurrence of these regions is an artifact 



5.2 The case with / = 1, 3, 5 

We now include the / = 5 contribution. The three equations 
obtained by projecting the GS equation (22) onto Z = 1, 3, 5 
are given in the Appendix A. The boundary conditions are 
essentially the same as in the previous Section; in particular, 
we look for the absolute minimum of asiRf + asiR)'^ , with 
fixed ai{R) = 1.93 • 10"=^ km. 

In Fig. 7 the profiles of the tetrad components of the 
magnetic field arc plotted for values of (q in the range < 
Co < 3.00 km^^. Fig. 8 shows the projections of the field lines 
in the meridional plane corresponding to the same values of 
Co- Comparing the results with the case Z = 1, 3 we see that 
the presence of the harmonic I = 5 modifies the magnetic 
field shape, but most of the features discussed in the previous 
Section are still present. 

An interesting difference is the following. While in the 
case Z = 1, 3 for ^ 0.40 km~^ we find field configurations 
which exhibit two disconnected regions where the function 
i/) has opposite sign and the magnetic field lines are confined 
(regions within dashed lines in Fig. 5), this does not occur 
when the I = 5 component is taken into account. This shows 
that the above feature has to be considered as an artifact of 
the truncation in the harmonic expansion. 



5.3 Higher order multipoles 

Up to now we have included components with Z < 7, neglect- 
ing the contribution from Z > 7. In order to test the accuracy 
of this approximation, we have studied the convergence of 
the harmonic expansion. To this purpose, we have solved 
the GS equation (22) including odd harmonic components 
up to Z = 7, for f = and Co = 0.61 km~^, and we have 
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Figure 9. The functions A'^^ (left panels) and A^^' (right panels) are shown for 
panels) in the meridional plane for < r < R. 



(upper panels) and Co = 0.61 km ^ (lower 




-1.0 -0.5 0.0 0.5 1.0 

x/R 

Figure 10. The projection of the field lines in the meridional plane is shown for Co = km~^ and Co = 0-30 km~^ respectively, and 
for I = 2,4. The dashed line corresponds to V = 0; the toroidal field is confined within the marked region. 



computed the quantities 

■'Pi<5{r,0) - il)i<3{r,e) 



ipi<7{r,9) - iJi<5{r,d) 



(36) 



These functions are shown in Fig. 9. They arc plotted only 
inside the star since outside they are much smaUer. Fig. 9 
shows that the error in neglecting I > 7, quantified by the 
function A^^^ is < 2% for Co = and < 4% for Co = 0.61 
km~^. Furthermore, a comparison of A(^) and A^'') shows 
that the harmonic expansion converges. 



5.4 An example of antisymmetric solution 

Here we show an example of a solution belonging to the 
antisymmetric family corresponding to I = 2,4. In Fig. 10 



we plot the field lines projected on the meridional plane, 

for Co = km^^ and Co = 0.30 km~^. Wo remark that 
the field lines arc antisymmetric with respect to the equa- 
torial plane; as a consequence, the total magnetic helicity 
is zero (see Section 6). Similar zero-helicity configurations 
have been considered in Braithwaite (2008a) . 



6 MAGNETIC HELICITY AND ENERGY 

The stationary configurations of magnetized neutron stars 
which we have found, depend on the value of the free pa- 
rameter Co, i.e. on the ratio between the toroidal and the 
poloidal components of the magnetic field. In this Section, 
we provide an argument to assign a value to Co. Furthermore, 
we compute the magnetic energy of the system to compare 
the contributions from poloidal and toroidal fields. 
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The total energy of the system (the star, the magnetic 
field and the gravitational field) can be determined by look- 
ing at the far field limit (r — > oo) of the spacetime metric 
(Misner et al. 1973; Thorne 1980). Following Colaiuda et al. 
(2008), loka & Sasaki (2004), we write the perturbed metric 
as 



ds^ = -e" [l + 2h{r, 6>)] dt^ + 



\ + —m{r,e) 
r 



dr 



[l + 2k{r, 6))] ^dO'^ + sinOdcj)^ 
+2't('r, e)dtdr + 2v{r, 6)dtd(j) + 2w{r, e)drd(t> (37) 

where, in particular, m{r,d) = ini{r) Pi {cos 9). The total 
mass-energy of the system is 



E = M + 5M 



(38) 



where M is the gravitational mass of the unperturbed star 
and 



5M = lim m()(r) . 



(39) 



In Appendix B, we discuss the equations which allow to 
determine E. We remark that 5M includes different contri- 
butions, due to magnetic energy, deformation energy, etc. 
In order to evaluate the magnetic contribution to E, 

it is convenient to use the Komar-Tolman formula (see for 
instance Straumann (2004), Chap. 4) for the total energy: 



(40) 



(where V is the 3-surface at constant time, is the time- 
like Killing vector, n** is the normalized, future directed nor- 
mal to V)] the magnetic contribution comes from the stress- 
energy tensor of the electromagnetic field, (3), i.e. 



'I 

Jv 

1 f°° 2 ^ r 

— r e dr 

2 Jo Jo 



mem \ ii v j-r r 

l^T g^^ ) ri'^n dV 



sine B^de. 



(41) 



The total (integrated) magnetic helicity Hm of the field 
configuration is 



Hm = 



(42) 



where is the f-component of the magnetic helicity 4- 
current, defined as 

1 



H° 



Explicitly, we have 

Hm = -2-K 



/ dr / [A, 
Jo Jo 



i),0 - iljAr,0\d9 



(43) 



(44) 



where 



e 2 



■V\o {W/i>\ - 1) • Q{]i>m - 1) 



Jo 



sin 6' 



Iv^/V-I - 1) 
■Q{\^/^\ - i)de' 



(45) 



The functional dependence of Hm on the potential of the 
toroidal field, Ar (see equation (44)), shows that regions of 



5M [10"" km] 
Em [10-^ km] 



0.5 



1.5 



2.5 



Figure 11. The functions SM and Em are plotted as functions 
of Co, for Z = 1, 3, 5 and Hm = 1.75 ■ IQ-^ W. 



space where the toroidal field vanishes do not contribute to 
the magnetic helicity. 

In ideal MHD, the magnetic helicity is a conserved 
quantity (Bekenstein 1987; Braithwaite & Spruit 2006). 
Thus, if we consider magnetic field configurations having the 
same value of the magnetic helicity and different energies, 
the lowest energy configuration is energetically favoured. 

In Fig. 11 we plot 6M and Em as functions of Co, for a 
fixed helicity Hm = 1.75 • 10~® km^. SM, and consequently 
the total energy, has a minimum at Co = 0.61 km~^. A fixed 
value of Hm corresponds, for any assigned value of ("o, to a 
different normalization constant B^oie- Since 5M, Hm and 
Em have the same quadratic dependence on the magnetic 
field normalization, this means that if we change Hm the 
plots of 5M and of Em as functions of Co are simply rescaled 
with respect to that shown in Fig. 11. Consequently, for any 
fixed value of Hm the position of the minimum of the total 
energy is the same as that shown in Fig. 11. We conclude 
that the configuration with Co — 0.61 km~^ is energetically 
favoured. This configuration is shown, among others, in Figs. 
7, 8. From Fig. 11 we also see that the contribution of the 
magnetic energy to 5M is of the order ~ 50%- 70%. 

In Fig. 12 we show the ratio of poloidal to total magnetic 
field energy, Ep/Em, as a function of Co, for the configura- 
tions (/ = 1,3,5) studied in this paper. This plot is inter- 
esting because the relative weight of the poloidal and the 
toroidal components of the field significantly affects many 
astrophysical processes involving magnetars, like magnetar 
activity (Woods & Thompson 2006), their thermal evolution 
(Pons et al. 2008), their gravitational wave emission (Cutler 
2002). It should be mentioned that the surface poloidal field 
is inferred from spin-down measurements which, however, 
provide no hint about the toroidal field hidden inside the 
star. We find that for Co = 0-61 knr\ Ep/Em ^ 0.93. 

In a recent paper (Braithwaite 2008b) the stability of 
magnetic field configurations of compact stars has been stud- 
ied in the context of Newtonian gravity, and assuming a 
polytropic equation of state. It has been found that axisym- 
mctric configurations are stable when 0.01 < Ep/Em ^ 0.8. 
Our configurations are outside this range, i.e. Ep/Em > 0.9. 
This difference may be attributed to several reasons: our 
models are computed in the framework of general relativ- 
ity, we use a more realistic equation of state, we choose a 
particular function C('0), which is linear in [\tp/'ip\ — l) (see 
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Figure 12. The ratio Ep/Em is shown as a function of Co, for 
I = 1,3,5. 



Eq. (10)). A different power-law dependence may lead to a 
diflferent contribution of the toroidal field. In any case, we 
observe a tendency in favour of models with predominant 
poloidal fields when using arguments of minimum energy, 
and wo do not think that other functional forms of may 
result in configurations with most of the energy stored in 
the toroidal field. This issue deserves further investigations. 



7 DISCUSSION AND CONCLUSIONS 



total energy at fixed magnetic hchcity. Wo find that, for our 
neutron star model (equation of state APR2, M — 1.4 Mq), 
this minimum occurs at Co = 0.61 km"\ Therefore, we ex- 
pect that in the early evolution of a strongly magnetized 
(fluid) neutron star, the natural final outcome after MHD 
equilibrium is estabhshod, arc twisted-torus configurations 
with geometries similar to our solutions. 

Finally we have computed the magnetic energy associ- 
ated to the poloidal and toroidal fields. We find that, al- 
though the amplitudes of both fields are of the same order 
of magnitude, and the toroidal field in the interior can be 
larger than the poloidal field at the surface (for instance, 
it is 2-3 times larger if Co = 0.61 km~^), the contribution 
of the toroidal field to the total magnetic energy is < 10%, 
because this field is non vanishing only in a finite region of 
the star. As mentioned in Section 6, a different power-law 
dependence of the function C(V') on [\ip/i>\ — l), may lead 
to a different contribution of the toroidal field and we plan 
to investigate this issue in a forthcoming paper. 
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In this paper wc find a twisted-torus family of solutions 
in the framework of general relativity. The toroidal compo- 
nent of the magnetic field vanishes outside the star: neither 
discontinuities (associated to surface currents), nor ad-hoc 
boundary conditions (total magnetic field vanishing outside 
the star) have been imposed; this is an improvement with 
respect to previous works (loka & Sasaki 2004; Colaiuda 
et al. 2008; Haskell et al. 2008). It should be stressed that 
twisted-torus configurations have boon found to emerge as a 
final outcome of Newtonian MHD simulations with generic 
initial conditions (Braithwaitc & Spruit 2004; Braithwaite 
& Nordlund 2006; Braithwaite & Spruit 2006). 

In order to have a twisted-torus configuration, there 
must be a non-linear relation between toroidal and poloidal 
fields, leading to couplings between different multipoles. We 
have investigated the contributions of different harmonics, 
and we have constructed equilibrium configurations with 
1 < i < 5. In order to fix the boundary conditions, we 
imposed that outside the star the dipolar component dom- 
inates, and minimized the Z > 1 contributions which, how- 
ever, remain non negligible. We find that there exist two 
particular, independent classes of solutions: those symmet- 
ric (with respect to the equatorial plane), in which all even 
order components vanish (02; = 0), and the antisymmet- 
ric solutions, characterized by the vanishing of odd compo- 
nents {a2i+i = 0). The latter have zero helicity by definition, 
therefore any solution minimizing energy at fixed helicity has 
a vanishing antisymmetric component. 

Our models also depend on a parameter, Co, which de- 
termines the ratio between toroidal and poloidal fields ( 
w Co^), and the length-scale of the region where the toroidal 
field is confined (oc C(7^)- Co increases, the amplitude of 
the toroidal field grows, but the region where it is confined 
shrinks. This parameter can be estimated by minimizing the 



APPENDIX A: GS EQUATION FOR THE CASE 
WITH Z = 1,3,5 

If we include the / = 1,3,5 components, the GS equation 
(22) projected onto the harmonics I = 1,3,5 gives the fol- 
lowing system: 

1 / -A »^ -x i^'-y , 2 



47r 



/"(3/4) Co (V' - SVIV-ZV^I + 2VVV5') 
Jo 



•e(|V'/'0| - l)sine de 



Co — -ci ai — -as 
5 \ 7 



{p + PV 



(Al) 



1 ( -X -x ^'-y , 12 



^ — l' PIT 



+^ j^^ (7/48)Co(v^-3V>IV>/^l+2VVV'') 

•e(|V'/V'l - 1)(3- 15cos^ 6')sin6'd6» 

= c,{p + Py{-a,--a^ + -aA , (A2) 



1 ( -X » -A t^'- A' , 30 

+^ j\ii/m) Co' (V' - 3V|V'/V5| + 2V'VV5') 

.9(1^/ Jl - i)(-315cos^^ + 210cos^^-15)^.^^^ 



8 



2 / 4 20 
■ ci{p + P)r [ ^as- — as 



(A3) 
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where 



-oi + 



03(3-15 008^0) 



+ 



a6(-315 cos* e + 210 cos^ 6 ■ 



15)" . 2, 
sm 



APPENDIX B: THE ENERGY OF THE 
SYSTEM 

The perturbation of the total energy of the system can be 
determined from the far field limit of the spacotime metric 
(Misner et al. 1973; Thorne 1980; loka & Sasaki 2004): 



5M = lim mo(r) , 



(Bl) 



where the perturbed metric is given by equation (37). The 
functions h{r,6) and m{r,8) are 



h{r,9) = ^/ii(r)Pi (cos 0) 
I 

m{r,e) = ^mi{r)Pi{cose) 



(B2) 



The perturbed Einstein equations {[tt] and [rr] compo- 
nents), projected onto I = 0, allow to determine the quantity 
mo(r): 

m'o - 'iTrr^-pjSpo = 

-(ffli) e +-(03)6 + r^ia^) e + 



72 ^2 , 450 2 
e~ 



^«3 + YI^a5 



/^\o(IV'/V'|-i)'e(|V'/V'|-i)^ 



h'o - e^^mo ( + 8nP ] - Anre^Spo = 



-{a[) +-{a',) +^(a^) - ^a, 
72e^ 2 450e^ 2 



7r3 
4r 



as 



llr3 



r Co (i^/^i 



i)^e(|v/V|-i)^d0 



(B3) 



where (5po is the / = component of the pressure 
perturbation (and vanishes outside the star), and = 
sin^^j^j 3 g aiPi 0. Using the relation (arising from T^Y^ = 
0) " ' 

6p'o = (^ + l)5p(.-(p + P)/io 



Co — -Cl I ffll — -03 



8 



10 



-^4(p + P)cjA„,__„3 + _a5 



10 , , „x / 4 20 
.-a5(p + P)ci(-a3--a5 



(B4) 



Eqns. (B3) can be rearranged in the form 



P' c 1/ /\2 -A , 6, /,2 -A , 15^ /^2 -A 

—5po = -{ai)e +-{a3)e + —{a^) e 

72 2 , 450 2 
+ 7737% 



/ 2 , 2 4J 

+5;:2«i + 7;:2«3 + YT 



Sp'o + 



^^\»(l'A/-0|-i)'e(|V'/^|-i)^o 

5po 



^(^ + l)+4nre\p + P) 



I 2A / 

-l-e mo( 



(P+-P)( 



4 / 3 ^ 
Co — -rCi I ai — -a3 



12 , / 2 8 10 \ 

10 , / 4 20 \ 1,2 6,2 
-YY«5Ci(^a3--a5l-^(ai) --{as) 



11 

15 

llr 

4r 



2e^ 2 72e^ 2 450e^ 2 
as 



iO , 2 , 26" 2 , 726-^ 2 , 



^\o' (IV'/V'I - i)'e(|V'/Vi| - i)^de 



By imposing a regular behaviour at r ~ we find 
mo{r^O)=Ar^ , 5po{r ^ Q) = Cr^ , 

where 

^ ^ 2(Pc + Pc)- (gf + aico) 
3 + 47r {r-^)^P. 



(B5) 
(B6) 



^ = 3 



(B7) 



The subscript c means that the quantity is evaluated at r — > 
0. We remark that the solution of (B5) does not depend on 
new arbitrary constants. Outside the star, the equation for 
mo reduces to 

° " 1 3 +7('*3) +^^(05) )'[}~~;rj 

2 2 , 72 2 ,450 2 
-^«i + ^«3 + 



mo 



450 2 
Ob 



(B8) 



+^«i + 7;3«3 + YI^ 
Solving (B5), (B8) we find 5M from (Bl) 
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